library(survey)
library(foreign)
library(ggplot2)
library(ggrepel)
library(cowplot)
library(dplyr)

## see wave1_cleaning.R, wave2_cleaning.R for data cleaning
load("~/Dropbox/Public_Conf_Mil/Data_and_code/wave2_data/clean_dataw2.RData")
load("~/Dropbox/Public_Conf_Mil/Data_and_code/wave1_data/clean_dataw1.RData")

## Create weighted survey design object
w1_design <- # wave 1 survey design object
  svydesign(
    id = ~ 1,
    weights = ~ weight,
    data = mil_conf.df
  )

w2_design <- # wave 2 survey design object
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = df
  )

addline_format <- function(x,...){
  gsub('\\s','\n',x)
}


#### TREATMENT EFFECTS (DIFF IN MEANS) ####

### Afghanistan, pandemic, riots

## Effective military - W1 - A6 - A6_d
## Ineffective military - W1 - A5 - A5_d
## Effective coronavirus - W2 - Assign1 = 3 - A1_3
## Ineffective coronavirus - W2 - Assign1 = 2 - A1_2
## NG positve - W2 - Assign1 = 4 - A1_4
## NG negative - W2 - Assign1 = 5 - A1_5

treat_df <- data.frame(conf = c(svyttest(Q8_d ~ ASSIGNMENT_A_6, w1_design)$estimate,
                                svyttest(Q8_d ~ ASSIGNMENT_A_5, w1_design)$estimate,
                                svyttest(Q11_d ~ A1_3, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_2, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_4, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_5, w2_design)$estimate),
                       ci_lo = c(svyttest(Q8_d ~ ASSIGNMENT_A_6, w1_design)$conf.int[1],
                                 svyttest(Q8_d ~ ASSIGNMENT_A_5, w1_design)$conf.int[1],
                                 svyttest(Q11_d ~ A1_3, w2_design)$conf.int[1],
                                 svyttest(Q11_d ~ A1_2, w2_design)$conf.int[1],
                                 svyttest(Q11_d ~ A1_4, w2_design)$conf.int[1],
                                 svyttest(Q11_d ~ A1_5, w2_design)$conf.int[1]),
                       ci_hi = c(svyttest(Q8_d ~ ASSIGNMENT_A_6, w1_design)$conf.int[2],
                                 svyttest(Q8_d ~ ASSIGNMENT_A_5, w1_design)$conf.int[2],
                                 svyttest(Q11_d ~ A1_3, w2_design)$conf.int[2],
                                 svyttest(Q11_d ~ A1_2, w2_design)$conf.int[2],
                                 svyttest(Q11_d ~ A1_4, w2_design)$conf.int[2],
                                 svyttest(Q11_d ~ A1_5, w2_design)$conf.int[2]),
                       cond = c(1:6))

fig5_1 <- ggplot(treat_df) + 
  geom_hline(yintercept = 0, colour = "red", lty = 1) + 
  geom_pointrange(aes(x = cond, y = conf, ymin = ci_lo, ymax = ci_hi), fatten = 2, size = 0.5) +
  geom_text(aes(x = cond, y = conf, label = sprintf("%0.2f", round(conf, digits = 2))), 
            size = 3,
            vjust = 2.5) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(1:6),
                     labels = addline_format(c("Won Afghanistan War", "Lost Afghanistan War",
                                               "COVID Success", "COVID Failure", 
                                               "NG Traditional", "NG Controversial"))) +
  ylab("Difference in Means from Control Condition") + ylim(-0.2, 0.2) +
  theme_bw()
fig5_1
ggsave("~/Dropbox/Public_Conf_Mil/New Chapter 5 Figures/effects_diff.jpeg", 
       plot = fig5_1, dpi = 700)

#### TREATMENT EFFECTS BY PARTY (DIFF IN MEANS) ####

treat2_df <- data.frame(rep = c(svyttest(Q8_d ~ A6_rep, w1_design)$estimate,
                                svyttest(Q8_d ~ A5_rep, w1_design)$estimate,
                                svyttest(Q11_d ~ A1_3_rep, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_2_rep, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_4_rep, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_5_rep, w2_design)$estimate),
                        dem = c(svyttest(Q8_d ~ A6_dem, w1_design)$estimate,
                                svyttest(Q8_d ~ A5_dem, w1_design)$estimate,
                                svyttest(Q11_d ~ A1_3_dem, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_2_dem, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_4_dem, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_5_dem, w2_design)$estimate),
                        ind = c(svyttest(Q8_d ~ A6_ind, w1_design)$estimate,
                                svyttest(Q8_d ~ A5_ind, w1_design)$estimate,
                                svyttest(Q11_d ~ A1_3_ind, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_2_ind, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_4_ind, w2_design)$estimate,
                                svyttest(Q11_d ~ A1_5_ind, w2_design)$estimate),
                        rep_lo = c(svyttest(Q8_d ~ A6_rep, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A5_rep, w1_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_3_rep, w2_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_2_rep, w2_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_4_rep, w2_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_5_rep, w2_design)$conf.int[1]),
                        rep_hi = c(svyttest(Q8_d ~ A6_rep, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A5_rep, w1_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_3_rep, w2_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_2_rep, w2_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_4_rep, w2_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_5_rep, w2_design)$conf.int[2]),
                        dem_lo = c(svyttest(Q8_d ~ A6_dem, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A5_dem, w1_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_3_dem, w2_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_2_dem, w2_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_4_dem, w2_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_5_dem, w2_design)$conf.int[1]),
                        dem_hi = c(svyttest(Q8_d ~ A6_dem, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A5_dem, w1_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_3_dem, w2_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_2_dem, w2_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_4_dem, w2_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_5_dem, w2_design)$conf.int[2]),
                        ind_lo = c(svyttest(Q8_d ~ A6_ind, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A5_ind, w1_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_3_ind, w2_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_2_ind, w2_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_4_ind, w2_design)$conf.int[1],
                                   svyttest(Q11_d ~ A1_5_ind, w2_design)$conf.int[1]),
                        ind_hi = c(svyttest(Q8_d ~ A6_ind, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A5_ind, w1_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_3_ind, w2_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_2_ind, w2_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_4_ind, w2_design)$conf.int[2],
                                   svyttest(Q11_d ~ A1_5_ind, w2_design)$conf.int[2]),
                        dcond = c(1,5,9,13,17,21),
                        icond = c(2,6,10,14,18,22),
                        rcond = c(3,7,11,15,19,23))

fig5_2 <- ggplot(treat2_df) + 
  geom_hline(yintercept = 0, colour = "black", lty = 1) +
  geom_pointrange(aes(x = rcond, y = rep, ymin = rep_lo, ymax = rep_hi), color = "red") +
  geom_pointrange(aes(x = icond, y = ind, ymin = ind_lo, ymax = ind_hi), color = "gray60") +
  geom_pointrange(aes(x = dcond, y = dem, ymin = dem_lo, ymax = dem_hi), color = "blue") +
  geom_text(aes(x = rcond, y = rep, label = sprintf("%0.2f", round(rep, digits = 2))), size = 3,
            vjust = 2) +
  geom_text(aes(x = icond, y = ind, label = sprintf("%0.2f", round(ind, digits = 2))), size = 3,
            vjust = 2) +
  geom_text(aes(x = dcond, y = dem, label = sprintf("%0.2f", round(dem, digits = 2))), size = 3,
            vjust = 2) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(2, 6, 10, 14, 18, 22),
                     labels = addline_format(c("Won Afghanistan War", "Lost Afghanistan War",
                                               "COVID Success", "COVID Failure", 
                                               "NG Traditional", "NG Controversial"))) +
  ylab("Difference in Means from Control Condition") +
  theme_bw()
fig5_2
ggsave("~/Dropbox/Public_Conf_Mil/New Chapter 5 Figures/effects_by_party.jpeg", 
       plot = fig5_2, dpi = 700)

#### ETHICAL TREATMENTS - WAVE 1 (DIFF IN MEANS) ####

trea3_df <- data.frame(conf = c(NA, svyttest(Q8_d ~ ASSIGNMENT_A_7, w1_design)$estimate,
                                svyttest(Q8_d ~ ASSIGNMENT_A_8, w1_design)$estimate, NA),
                       ci_lo = c(NA, svyttest(Q8_d ~ ASSIGNMENT_A_7, w1_design)$conf.int[1],
                                 svyttest(Q8_d ~ ASSIGNMENT_A_8, w1_design)$conf.int[1], NA),
                       ci_hi = c(NA, svyttest(Q8_d ~ ASSIGNMENT_A_7, w1_design)$conf.int[2],
                                 svyttest(Q8_d ~ ASSIGNMENT_A_8, w1_design)$conf.int[2], NA),
                       cond = c(0.5,1,2,2.5))

fig5_3 <- ggplot(trea3_df) + 
  geom_hline(yintercept = 0, colour = "red", lty = 1) + 
  geom_pointrange(aes(x = cond, y = conf, ymin = ci_lo, ymax = ci_hi), fatten = 2, size = 0.5) +
  geom_text(aes(x = cond, y = conf, label = sprintf("%0.2f", round(conf, digits = 2))), 
            size = 3,
            vjust = 2.5) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(1,2),
                     labels = addline_format(c("Ethical Generals", "Unethical Generals"))) +
  ylab("Difference in Means from Control Condition") + ylim(-0.3, 0.1) +
  theme_bw()
fig5_3
ggsave("~/Dropbox/Public_Conf_Mil/New Chapter 5 Figures/ethical_diff.jpeg", 
       plot = fig5_3, dpi = 700)

#### ETHICAL TREATMENTS BY PARTY (DIFF IN MEANS) ####

treat4_df <- data.frame(rep = c(svyttest(Q8_d ~ A7_rep, w1_design)$estimate,
                                svyttest(Q8_d ~ A8_rep, w1_design)$estimate),
                        dem = c(svyttest(Q8_d ~ A7_dem, w1_design)$estimate,
                                svyttest(Q8_d ~ A8_dem, w1_design)$estimate),
                        ind = c(svyttest(Q8_d ~ A7_ind, w1_design)$estimate,
                                svyttest(Q8_d ~ A8_ind, w1_design)$estimate),
                        rep_lo = c(svyttest(Q8_d ~ A7_rep, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A8_rep, w1_design)$conf.int[1]),
                        rep_hi = c(svyttest(Q8_d ~ A7_rep, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A8_rep, w1_design)$conf.int[2]),
                        dem_lo = c(svyttest(Q8_d ~ A7_dem, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A8_dem, w1_design)$conf.int[1]),
                        dem_hi = c(svyttest(Q8_d ~ A7_dem, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A8_dem, w1_design)$conf.int[2]),
                        ind_lo = c(svyttest(Q8_d ~ A7_ind, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A8_ind, w1_design)$conf.int[1]),
                        ind_hi = c(svyttest(Q8_d ~ A7_ind, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A8_ind, w1_design)$conf.int[2]),
                        dcond = c(1,5),
                        icond = c(2,6),
                        rcond = c(3,7))

fig5_4 <- ggplot(treat4_df) + 
  geom_hline(yintercept = 0, colour = "black", lty = 1) +
  geom_pointrange(aes(x = rcond, y = rep, ymin = rep_lo, ymax = rep_hi), color = "red") +
  geom_pointrange(aes(x = icond, y = ind, ymin = ind_lo, ymax = ind_hi), color = "gray60") +
  geom_pointrange(aes(x = dcond, y = dem, ymin = dem_lo, ymax = dem_hi), color = "blue") +
  geom_text(aes(x = rcond, y = rep, label = sprintf("%0.2f", round(rep, digits = 2))), size = 3,
            vjust = 2) +
  geom_text(aes(x = icond, y = ind, label = sprintf("%0.2f", round(ind, digits = 2))), size = 3,
            vjust = 2) +
  geom_text(aes(x = dcond, y = dem, label = sprintf("%0.2f", round(dem, digits = 2))), size = 3,
            vjust = 2) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(2, 6),
                     labels = addline_format(c("Ethical Generals", "Unethical Generals"))) +
  ylab("Difference in Means from Control Condition") +
  theme_bw()
fig5_4
ggsave("~/Dropbox/Public_Conf_Mil/New Chapter 5 Figures/ethical_by_party.jpeg", 
       plot = fig5_4, dpi = 700)
